Cell Systems

Predicting single-cell transcriptional responses to perturbations is central to building the virtual cell, yet recent benchmarks show that simple baseline methods often outperform complex models, and model comparisons depend on the evaluation metric. Most studies assume that preprocessed RNA sequencing data are reliable ground truth for both training and evaluation. Here, we test this assumption by measuring the reliability of perturbations and their alignment with shared perturbation responses, classifying each perturbation as specific, shared, or unreliable. Among 7,170 perturbations from 29 datasets, 65% are unreliable, 11% shared, and 24% specific. Applying these quality labels to published benchmarks shows that model comparisons depend on perturbation quality. Training with reliable perturbations alone matches or outperforms full-data performance while using 55% of all training perturbations. Our framework also enables prospective experimental design: for most perturbations, a 28-cell pilot experiment accurately predicts how many cells a full screen needs to be reliable.

Predictive logical models have made substantial contributions to biomedical research by facilitating the identification of regulatory mechanisms driving biological processes. Yet, their construction remains a major bottleneck, as it typically relies on extensive prior knowledge at the genome scale. The emergence of high-throughput omics technologies has opened new opportunities for data-driven logical model construction by providing high-resolution observations of biological states. In parallel, formal constraint-based approaches have emerged, reformulating executable model inference as a constraint-satisfaction problem. However, transforming omics observations into biologically meaningful Boolean state abstractions and dynamical constraints suitable for logical synthesis engines remains challenging, particularly for poorly understood biological systems and multi-condition experimental settings. To address these challenges, we introduce scBOLT, a semi-automated framework for Boolean constraint engineering and logical model inference from multi-condition single-cell transcriptomes. Beyond these methodological contributions, scBOLT substantially improves computational tractability, enabling the inference of Boolean networks comprising several hundred components within practical runtimes.

Genetic perturbations can reshape cell populations by altering the relative abundance of specific cell types and states within the profiled population, including increases, decreases and states that become detectable after perturbation. Pooled single-cell screens, such as Perturb-seq, measure such responses at scale. However, only a small fraction of possible perturbations can be tested experimentally. A central challenge is therefore to predict compositional shifts induced by unseen perturbations. Many perturbation-prediction methods do not directly optimize for this outcome; instead, they predict gene-expression responses and infer cell-type and cell-state composition downstream. Surprisingly, we find that even models that accurately predict perturbation-induced changes in average gene expression perform poorly at forecasting these compositional shifts. To address this gap, we present Pop-Corn, a method that directly predicts how a perturbation reshapes cell-type composition without reconstructing gene expression. In the primary T-cell benchmark, Pop-Corn predicted the overall cell-state composition of held-out perturbations more accurately than the evaluated expression-prediction pipelines, while better preserving the diversity of observed cell states. We further extend Pop-Corn to intact tissue, where it predicts perturbation-induced cell-type proportion changes in local cellular neighborhoods and uses attention patterns to generate hypotheses about context-dependent cellular interactions. Retrospective virtual screens support the use of Pop-Corn to prioritize perturbations for experimental follow-up according to their predicted effects on cell-state composition.

Machine learning models in single-cell biology increasingly forecast differentiation, reprogramming and therapeutic response from early transcriptomic profiles. Testing whether a model has learned real biology requires held-out cells. Single-cell data, however, are grouped: cells from the same clone, patient or batch share the same label. A random split therefore places relatives of each test cell, carrying its label, in the training set, and a model can score well by memorizing a relative instead of learning a transferable rule. Grouped validation removes this leakage but leaves far fewer independent units behind each error bar. Here we show how to estimate this leakage before training any model, from two properties of the data: exposure, the fraction of test cells with relatives in training, and retrievability, how often a nearest-neighbor search returns such a relative rather than an unrelated cell. Across lineage-barcoded and patient data, exposure determines whether a random split opens a leakage channel, and retrievability determines how much it can inflate the score. The inflation is negligible where cell state has decoupled from ancestry, much larger where clonal sisters remain close in expression space, and in a patient cohort large enough to overturn a clinical conclusion. We also provide eakcheck, which computes both properties in seconds, before the outcome model is fitted.

It is widely accepted in genomics that coexpression of RNA transcripts suggests a commonality of function. This intuition is explicitly leveraged in machine learning methods that predict gene function, where it is often combined with other features such as protein interactions and sequence similarity. For example, including coexpression data from human tissue expression boosts performance for predicting Gene Ontology annotations. However, the biological underpinnings of this observation have not been well-investigated. Building on earlier results from our group, in this work we show that gene function is predictable from coexpression substantially because it reflects differences in expression between cell types, and these differences are also intrinsic to the ground truth labels. Using simulations and analyses of real data, we show that variance in the cellular composition of bulk samples impacts function learnability and attribute this to cell type marker gene content in the GO terms. We further show that cell type profiles, where the relationship between gene expression and cell type is made transparent, are effective for predicting gene function while increasing interpretability. These results indicate that function prediction models trained on bulk coexpression are largely limited to cell-type-level resolution rather than fine-grained biochemical function, with direct consequences for how such predictions should be interpreted.

Predicting how a genetic perturbation reshapes a cells transcriptome is a central goal of computational biology. Previous studies report that the mean response across training perturbations rivals specialized models on standard accuracy metrics, even though it cannot distinguish which perturbation occurred. Across 2,270 CRISPRi perturbations measured in each of six cell lines, we show that this apparent paradox reflects a conserved organization of perturbation responses. Perturbations span a continuum from responses strongly aligned with the mean to more targeted responses that depart from it. Crucially, a perturbations position along this continuum is conserved across cell lines (Kendalls W = 0.59) and predictable from STRING protein-interaction embeddings (R2 = 0.35). We formalize this structure with COMPASS, an interpretable linear model that decomposes each response into shared and gene-specific components and estimates them separately. The shared-response component is modeled as a cell-line-wide response scaled by a perturbation-specific coefficient. This coefficient is strongly conserved across cell lines. The residual gene-specific component--which is moderately conserved across cell lines--recovers pathway-level programs. COMPASS outperforms scGPT, CPA, GEARS, GenePert, and SO_SCPLOWTATEC_SCPLOW in both response accuracy (de-biased Pearson delta 0.34 vs. [≤] 0.32) and perturbation discrimination (cosine PDS gain 0.23 vs. [≤] 0.08). These results recast perturbation prediction across cellular contexts as component-wise inference, with each component estimated from the evidence best suited to it.

Combination therapies are central to cancer treatment, but exhaustive screening is impractical. We introduce PHAROS, a framework that turns a pretrained single-cell perturbation model into a target-directed search engine for drug combinations. PHAROS predicts how a cell population changes under a drug, one drug at a time, then chains these predictions together to simulate drug combinations. It scores each simulated outcome against the desired target state and uses a search algorithm to find the most promising combinations, all without retraining the underlying model. Across two independent combinatorial perturbation datasets, PHAROS recovered exact or mechanism-matched two-drug responses in cell lines, both seen and unseen during model training. Its rankings were specific to the requested conversion and were not explained by single-drug effects, additive effects, or shared mechanism of action. In exploratory analyses of patient-derived metastatic HR+/HER2$-$ breast tumors and basal cell carcinoma (BCC), PHAROS prioritized FDA-approved regimens, distinguished combinations by their predicted tumor-versus-immune objective profiles, and nominated pathway-level hypotheses, while explicitly identifying both tumor cohorts as outside the model's supported distribution. PHAROS provides a modular route from pretrained virtual-cell models to inverse, single-cell combination-screening platforms.

Mechanistic models of biological networks predict but cannot be calibrated at scale, while learned representations describe but expose no dynamical operator. Both lack an inverse step mapping data onto the few dynamical quantities the data determine. We develop Unified Network Spectroscopy (UNS) as that operator and establish when it is well posed. Reading fluctuations about a physiological setpoint through the Lyapunov relation of the linear noise approximation, UNS returns dominant timescales and frequencies, the subspace on which they act, and the projected noise. The inversion itself is established; our contribution is architectural. Three results govern its success. First, biochemical noise is generically anisotropic, and extrinsic modulation makes the dynamics non-normal. Assuming isotropy can invert mode rankings and destroy precision-matrix topology recovery; we give the exact, checkable condition under which recovery survives. Second, equal-time covariance identifies neither the dynamics nor their spectrum: the residual freedom is precisely the probability current that non-equilibrium living systems generically carry. Third, population snapshots replace time series only under ergodicity; heterogeneity biases inferred relaxation rates downward, making systems appear closer to bifurcation, whatever the sample size. Lagged covariance or per-cell velocity restores identifiability and turns noise structure into a recovered, falsifiable quantity. We provide a three-regime ergodicity taxonomy with diagnostics and show that the variational operator remains well posed when samples are fewer than variables. The output is a typed interface contract for hierarchical model composition, carrying its uncertainty, attainment and validity domain, with controllability inheritance as a falsification test.

Microbiome measurements are commonly judged by how faithfully they reconstruct an underlying biological profile, using taxonomic similarity, prediction accuracy, or compositional distance. Yet faithfulness and usefulness are not the same: a distorted measurement may still support the correct intervention, while a small error can be decisive when it falls near the boundary between two treatment choices. We therefore evaluate measurement quality by whether information loss changes the decisions and outcomes that remain attainable. We treat the intervention-relevant state of the gut community as a hidden quantity seen only through an imperfect, information losing measurement, mapped by a decision rule to a graded intervention. Each decision rule is evaluated by net host benefit and ecological margin, generating a feasible set of attainable outcome pairs. We define the measurement-induced decision loss (MIDL) as the maximum net host benefit lost under imperfect measurement relative to perfect observation. For scalar affine information loss, increasing degradation is a Blackwell garbling; under post-processing closure, more informative measurement weakly expands the feasible set and cannot reduce optimal host benefit. In the binary-state model, MIDL increases linearly while an observation-dependent decision rule remains optimal and plateaus once the optimal policy becomes signal-independent. We extend the framework to finite states, state-specific loss, multi-segment stool aggregation, and ecological-margin constraints. The framework distinguishes measurement faithfulness from decision usefulness and provides a principled basis for comparing stool and richer measurements by the interventions and trade-offs they support.

Pooled Perturb-seq screens are often interpreted through a linear steady-state model in which each targeted gene supplies an additive perturbation input. Although this formulation identifies an operator on the observed response subspace in closed form, identifiability alone does not establish that inferred operators generalize to new perturbation targets. We evaluated fixed program-space additive inputs in three Perturb-seq screens -- Replogle K562 essential, Replogle RPE1 essential, and Jost 2020 -- using target-grouped nested cross-validation and matched signal-to-noise linear-truth controls. For inverse operator prediction, real held-out error remained at the predict-zero baseline (0.96, 1.00, and 1.00), whereas matched linear controls achieved 0.18, 0.53, and 0.72. This discrepancy persisted across random K562 target panels, an overdetermined Jost dimension sweep, direction-only fits, ridge and truncated-SVD estimators, and five linear ground-truth ensembles. Forward target-held-out prediction was more nuanced: fixed, footprint, and learned linear encodings modestly beat a predict-training-mean baseline in K562 but not in RPE1 or Jost, and the K562 gains recovered only about 7% of the gap to the matched linear control. On the tested screens, fitted fixed program-space operators should not be interpreted as validated target-generalizing regulatory operators. We release anchor-op, a reproducible toolkit for target-grouped evaluation, identifiability checks, and matched-geometry positive controls for perturbation-response models.

Inductive signals can direct cell fate decisions, yet the diversity of cell identities far exceeds the available signaling repertoire. While developing systems resolve this paradox by applying temporal sequences of minimalistic signals to hierarchically organized gene regulatory networks (GRNs), the rules governing robust sequential fate navigation remain largely unknown. To address this, we modeled how sequential signals steer trajectories toward arbitrary fates within a Waddington-like manifold, enumerating parameters and combinatorial logics of cross-inhibition and self-activation (CIS) modules, and various ways of applying the sequential signals. We uncover a fundamental conflict between commitment stability and inductive plasticity that severely limits sequential fate navigation. Crucially, introducing non-inductive gap intervals systematically resolves this bottleneck by dissipating kinetic leakage without sacrificing upstream stability. Finally, we demonstrate that polarization-division cycle can autonomously generate signals meeting these temporal requirements. Together, our findings provide a minimal model demonstrating the feasibility and design principles for driving hierarchical cell fate decisions through sequential pulses of a restricted signaling repertoire.

Virtual cell models aim to predict how cells respond to perturbations, yet accurate prediction of post-perturbation gene expression can conceal failures to recover context-specific responses. Here we identify context geometry compression (CGC), in which models preserve cellular identity and shared perturbation effects while collapsing differences in how the same perturbation acts across contexts. Across biological systems, reproducible context-specific responses remained difficult to recover; in a large drug-response matrix, this persisted even with nearly complete perturbation coverage. The same drug-response predictions achieved higher normalized recovery in predefined pathway summaries than in gene-level response profiles, while coarse response programs were more predictable in human lymphoblastoid cell lines. Compact baseline transcriptomic representations retained most full-RNA predictive performance for patient-specific ex vivo drug responses. Together, these findings establish recoverable biological resolution, rather than measurement dimensionality, as a criterion for evaluating virtual cell predictions and matching model outputs to the biological detail supported by evidence.

Identifying which gene drives a desired cell-state transition is a foundational inverse problem in target discovery and cell-state engineering, yet far less developed than forward prediction. We show that prevailing benchmarks inflate performance through two artefacts: training leakage, where test-perturbation cells are seen during fitting, and evaluation leakage, where the perturbed gene's own expression change exposes the answer. Here, we introduce Cause-excluded Effect Evaluation (CEE), a framework that separates train and test perturbation targets, removes target-gene self-effects from expression profiles and scores candidates in downstream response space. With CEE, we find that simple structure-based methods consistently outperform complex deep learning models across 16 single-cell perturbation datasets, revealing that the recoverable signal is governed by low-rank and additive structure in perturbation responses rather than by model complexity. Our work disentangles evaluation artefacts from genuine predictive capacity and offers a design principle for developing perturbation target identification methods.

During development, cells traverse gene expression states where their local dynamical sensitivity changes sharply, yet existing computational methods provide limited access to when and where this sensitivity peaks along a trajectory. Here we introduce scJDO (single-cell Jacobian Differential Operators), a framework that characterizes how local dynamical sensitivity evolves during cell fate transitions, together with an explicit account of what that representation can and cannot recover from snapshot data. scJDO treats time-indexed Jacobians as explicit analytical objects, projecting the temporal sequence of operators into a shared subspace and decomposing it into recurrent operator archetypes with interpretable temporal activation profiles. Unlike methods that derive Jacobians from splicing-kinetic vector fields, scJDO learns a neural drift field directly from cell-state geometry via diffusion score matching, enabling Jacobian analysis on trajectory-resolved scRNA-seq datasets regardless of splicing-data availability. Applied to a dense time-course of induced pluripotent stem cell (iPSC) reprogramming, scJDO resolves a quantitative operator-level signature that distinguishes diverted from productive fate: the productive trajectory executes a sequential handoff from an early MEF-exit operator regime to a late pluripotency-associated regime, whereas the diverted trajectory maintains the early regime and instead activates a distinct stress-associated archetype. We validate scJDO across four settings: synthetic benchmarks with analytically known ground truth, branching hematopoiesis, dense real time-course reprogramming, and a perturbational setting using Schrodinger bridges in K562 CRISPRi Perturb-seq. We compare against the two most widely used single-cell Jacobian methods on a dataset where all three are runnable, finding that scJDO shares significantly more gene-level and directional operator structure with Dynamo than expected by chance while providing operator-level analysis on datasets without splicing kinetics. We further characterize the boundary of the representation directly. At a fate-decision saddle, eight mathematically distinct readouts of the same learned drift field are consistent with a single explanation: a drift field fit to snapshot density reproduces density-dominant separation between committed branches rather than the low-variance transverse instability that defines the decision. Together, scJDO provides an operator-level view of single-cell dynamics and an explicit characterization of its own identifiability boundary.

Predictive models of cellular perturbation are often judged by how closely they reconstruct molecular states after unseen interventions. We show that high state-level similarity can coexist with loss of the relationships that distinguish perturbations, a failure we term Intervention Geometry Compression (IGC). Across established models and perturbation settings, unseen interventions show weakened global and local geometry, reduced between-intervention variance and spectral collapse. The failure is not primarily explained by response-space capacity. Instead, diagnostic projections localize much of the missing geometry to a small number of residual response directions learned from seen interventions; these directions outperform complexity-matched random subspaces and replicate in an independent Jiang perturbation resource. Polarity captures part, but not all, of this continuous orientation signal. Time-resolved analyses further show that correct trajectory entry markedly improves downstream propagation, while a held target's own early empirical response rapidly reveals endpoint orientation. Finally, same-target empirical anchoring transfers intervention identity across contexts far more effectively than increasing exposure to other interventions. These results identify intervention-coordinate assignment as an information bottleneck in virtual-cell generalization and support a design principle: empirically anchor intervention identity, then use models to generalize anchored effects across cellular contexts.

Predicting single-cell responses to genetic perturbations could reveal the vast combinatorial space of perturbations and cellular contexts that is infeasible to measure experimentally, yet current deep learning models generalize poorly and often fail to outperform simple baselines. Here we demonstrate that generalizability in perturbation prediction requires identifying and representing distinct components of cellular response rather than on increasing model complexity alone. We introduce a decomposition framework that explicitly separates transcriptional responses into global, perturbation-specific, cell-line-specific, and perturbation-by-cell-line interaction components. Applied to four CRISPR-interference Perturb-seq screens on multiple cell lines, our framework reveals that these components have distinct structures and information requirements. The global response component is low-dimensional, reflects recurrent proliferation and stress response programs, and can be inferred from control gene expression. In contrast, the perturbation and cell-line specific components are high-dimensional and cannot be recovered from control expression. We therefore develop response-component-aligned models that map biological priors, such as gene coessentiality, onto the geometry of observed transcriptional responses. Critically, this alignment enables simple linear or multilayer perceptron (MLP) based models to outperform state-of-the-art architectures across multiple generalization settings, including unseen cell lines and combinations of unseen perturbations. Together, our framework for response decomposition and alignment provides a principled basis for evaluating and designing perturbation-prediction models, showing that generalization depends primarily on matching biological information to the response components rather than on model complexity alone.

Quantitatively relating transcription factor (TF) input to gene expression output is central to understanding mammalian gene regulation and essential for designing predictable synthetic expression systems. However, even minimal synthetic systems often exhibit unexplained behaviors. In a widely used inducible mammalian expression system, we show that transcriptional responses appear graded and sigmoidal at the population level but are largely all-or-none at the single-cell level. By combining single-cell sorting and single-molecule footprinting with mathematical modeling of transcriptional regulation, we found that this behavior is not caused by bursty transcription or bistability, but by long-lived, chromatin-encoded variability in TF occupancy and activation strength. This variability produced a range of activation thresholds in switch-like single-cell responses that were stable over time, resulting in bimodal gene expression across the population. These results advance our basic understanding of how TFs interact with chromatin to modulate quantitative features of single-cell and population level transcriptional responses.

In high-throughput single-cell transcriptomics (p {approx} 20,000 genes), performing double machine learning (DML) causal inference on q {approx} 5,000 target genes requires nuisance function fits that grow linearly with the number of targets (Kf cross-fitting folds, Kf = 5 or 10), far exceeding feasible computational budgets, especially with deep learning. We propose a Randomized Partition Strategy (RPS): randomly divide target genes into groups, share one background compression per group, reducing deep learning model training to q/m runs (m = group size) --- a factor of m savings. The cost of grouping is accuracy loss --- we prove that the cumulative deviation of the estimator follows a one-dimensional drift-free symmetric random walk, with diffusion variance growing linearly with group size and mean squared displacement equaling the mean squared error, so accuracy loss is predictable: m = 1 is always optimal, accuracy cost is monotonically increasing, and a small accuracy sacrifice yields m-fold compute savings. On GSE189050 SLE single-cell data (Memory B cells, n = 2120), both PCA and DL methods converge to the same conclusion, confirming the random walk mechanism is method-independent; an unexpected finding is that DL diffusion growth is only 16%, far slower than PCAs 7.4 times. This work provides a quantifiable theoretical foundation for compute strategy selection in single-cell high-dimensional causal inference.

Single-cell perturbation experiments measure how interventions alter cellular phenotypes. However, the number of possible perturbations and biological contexts far exceeds what can be tested experimentally. Motivated by this constraint, predictive models aim to extrapolate cellular responses to unseen conditions. Despite substantial efforts in model development, benchmark studies have reached inconsistent conclusions about the capabilities of current perturbation-response models. A major challenge is that evaluation protocols vary widely across studies, making results difficult to compare. Furthermore, the lack of consensus on evaluation hampers progress because it is unclear which predictive capabilities new models should prioritize. To help build consensus on evaluation principles, we develop a taxonomy that decomposes evaluation protocols into their representation, metric, score transformation, and reporting strategies. We characterize how these choices determine which aspects of prediction quality a benchmark measures and discuss criteria for selecting and assessing protocols in relation to specific benchmarking goals. We additionally provide scPertEval, a Python package with reference implementations of selected evaluation protocols, and use it to assess protocol behavior across seven publicly available single-cell perturbation datasets. By making evaluation choices and their underlying trade-offs explicit, we aim to stimulate a community discussion about developing more comparable and task-aligned evaluation protocols. Graphical Abstract O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=129 SRC="FIGDIR/small/740433v1_ufig1.gif" ALT="Figure 1"> View larger version (32K): org.highwire.dtl.DTLVardef@bd1e3corg.highwire.dtl.DTLVardef@c26c2org.highwire.dtl.DTLVardef@1c49afdorg.highwire.dtl.DTLVardef@9b988c_HPS_FORMAT_FIGEXP M_FIG C_FIG

Biological AI learns primarily from existing observations, but many questions cannot be answered from available data alone. Here we show that AI can instead acquire knowledge by acting directly on biological systems and learning from the consequences. We developed a closed-loop framework in which autonomous agents design protein variants, construct and characterize them in a robotic laboratory, learn from the resulting experimental feedback, and decide what experiments to perform next. We then allowed the system to operate continuously and without human intervention for approximately one month, during which multiple agents independently explored protein sequence space while learning from shared experimental experience. Applied to glycoside hydrolases, the agents discovered enzymes with substantially altered substrate specificity toward non-native sugars and progressively learned the structure of the underlying sequence-function landscape. The resulting experimental experience also revealed determinants of substrate specificity and protein expression that were not specified as learning objectives. These results demonstrate that AI can autonomously interact with biology over extended periods to acquire knowledge through experience, establishing a framework for biological discovery driven by continuous experimental interaction.